function [Ex, rho_a, phase] = MT1D_FDM(Length, Nz, S, fre)
% MT1D_FDM 一维大地电磁正演模拟（有限差分法）
% 用于计算层状介质的大地电磁响应
%
% 输入参数：
%   Length  - 计算区域的深度（单位：米）
%   Nz      - 网格剖分单元数
%   S       - 电导率分布向量（单位：S/m），长度应为Nz+1
%
% 输出参数：
%   Ex      - 电场响应矩阵（Nz+1 × 频率数）
%   rho_a   - 视电阻率向量（单位：Ω·m）
%   phase   - 相位向量（单位：度）
%
% 参考文献：
%   偏微分方程的有限差分法及地球物理应用，第 6章

% 物理常数定义
mu = 4e-7 * pi;        % 真空磁导率 (H/m)
eps = 8.8419e-012;     % 真空介电常数 (F/m)，代码中未使用但保留

% 网格参数计算
Dz = Length / Nz;      % 网格间距

% 预分配输出变量
Ex = zeros(Nz+1, length(fre));    % 电场响应矩阵
rho_a = zeros(1, length(fre));    % 视电阻率
phase = zeros(1, length(fre));    % 相位

% 频率循环：对每个频率点独立计算
for i = 1:length(fre)
    omega = 2 * pi * fre(i);      % 角频率
    
    % 初始化稀疏矩阵和右端项
    P = sparse(Nz+1, 1);          % 右端项向量（源项）
    K = sparse(Nz+1, Nz+1);       % 系数矩阵（刚度矩阵）
    
    % 构建内部节点的有限差分方程
    % 方程形式：d²Ex/dz² + iωμσEx = 0
    for j = 2:Nz
        % 中心差分格式，考虑电导率在节点间的平均值
        sigma_avg = (S(j-1) + S(j)) / 2;  % 界面电导率平均
        
        % 设置三对角矩阵元素
        % 左邻居：1/Dz²
        % 自身：iωμσ - 2/Dz²  
        % 右邻居：1/Dz²
        K(j, j-1:j+1) = [1/Dz^2, ...
                         sqrt(-1)*omega*mu*sigma_avg - 2/Dz^2, ...
                         1/Dz^2];
    end
    
    % 上边界条件（地表）：Dirichlet边界条件
    % Ex(0) = 1 (归一化电场)
    K(1, 1) = 1;    % 单位矩阵元素
    P(1) = 1;       % 单位源项
    
    % 下边界条件（底部）：Robin边界条件
    % dEx/dz + a*Ex = 0，其中a为衰减系数
    a = sqrt(-sqrt(-1)*omega*mu*S(end));  % 波数（衰减系数）
    K(end, end-1) = -1/Dz;        % 一阶导数项
    K(end, end) = 1/Dz + a;       % 边界条件主项
    P(end) = 0;                   % 无源项
    
    % 求解线性方程组：K * Ex = P
    Ex(:, i) = K \ P;             % 直接法求解
    
    % 计算地表电磁场分量用于视电阻率和相位计算
    Ex_g = Ex(1, i);              % 地表电场
    
    % 计算地表磁场（Hy）：使用四点差分公式提高精度
    % 公式：( -11*Ex1 + 18*Ex2 - 9*Ex3 + 2*Ex4 ) / (6*Dz)
    Hy_g = (-11*Ex(1, i) + 18*Ex(2, i) - 9*Ex(3, i) + ...
            2*Ex(4, i)) / (6 * Dz) / (sqrt(-1) * mu * omega);
    
    % 计算波阻抗 Z = Ex / Hy
    Z = Ex_g / Hy_g;
    
    % 计算卡尼亚视电阻率：ρ_a = |Z|² / (ωμ)
    rho_a(i) = abs(Z)^2 / (mu * omega);
    
    % 计算阻抗相位（转换为角度制）
    phase(i) = -atan2(imag(Z), real(Z)) * 180 / pi;
end

end